Perturbation theory vs. simulation for tadpole improvement factors in pure gauge 
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We calculate the mean link in Landau gauge for Wilson and improved SU(?>) anisotropic gauge 
actions, using two loop perturbation theory and Monte Carlo simulation employing an accelerated 
Langevin algorithm. Twisted boundary conditions are employed, with a twist in all four lattice 
directions considerably improving the (Fourier accelerated) convergence to an improved lattice Lan- 
dau gauge. Two loop perturbation theory is seen to predict the mean link extremely well even into 
the region of commonly simulated gauge couplings and so can be used remove the need for numer- 
ical tuning of self-consistent tadpole improvement factors. A three loop perturbative coefficient is 
inferred from the simulations and is found to be small. We show that finite size effects are small 
and argue likewise for (lattice) Gribov copies and double Dirac sheets. 

PACS numbers: 11.15.Ha,12.38.Gc 



I. INTRODUCTION 

Tadpole improvement is now widespread in lattice field 
theory [1]. Without it, lattice perturbation theory begins 
to fail on distance scales of order 1/20 fm. Perturba- 
tion theory in other regularisations, however, seems to 
be phenomenologically successful down to energy scales 
of the order of 1 GeV (corresponding to lattice spacings 
of 0.6 fm) [2]. 

The reason is that the bare lattice coupling is too small 
[1, 2]. To describe quantities dominated by momenta of 
order the cut-off scale (tt/o), it is appropriate to expand 
in the running coupling, as, evaluated at that scale. The 
bare coupling, however, deviates markedly from this and 
its anomalously small value at finite lattice spacing can 
be associated with tadpole corrections [2] . These tadpole 
corrections are generally process independent and can be 
(largely) removed from all quantities by modifying the 
action. This corresponds to a resumming of the pertur- 
bative series to yield an expansion in powers of a new, 
"boosted" coupling that is much closer to as{-K/a). 

Perturbatively this amounts to adding a series of radia- 
tive counterterms to the action. Such a series is obtained 
by dividing each gauge link in the action by an appro- 
priate expansion, u'^^"^\ It is sufficient that this series is 
known only up to the loop order of the other quantities 
we wish to calculate using the action. 

The factor m*^^"""' is not unique, but it should clearly be 
dominated by ultraviolet fluctuations on the scale of the 
cut-off. The two most common definitions are the fourth 
root of the mean plaquette (a gauge invariant definition) 
and the expectation value of the link in Landau gauge. 
Both are successful, although lattice results suggest some 
arguments for preferring the latter [3]. In this paper we 
discuss Landau gauge mean link tadpole improvement. 

For Monte Carlo simulations, each gauge link in the ac- 
tion is divided by a numerical factor, u^^'~^\ Its value is 
fixed by a self-consistency criterion; the value measured 
in the simulation should agree with the parameter in 



the action. Obtaining such numerical agreement requires 
computationally expensive tuning of the action. In many 
cases this cost is prohibitive, such as the zero tempera- 
ture tuning of highly anisotropic actions for use in finite 
temperature simulations. As non-perturbative phenom- 
ena should not affect the cut-off scale, a sufficiently high 
order perturbative series should predict m*^'^^^ such that 
the subsequent numerical tuning is unnecessary. 

In this paper we present the tadpole improvement fac- 
tors calculated to two loop order using lattice perturba- 
tion theory. This covers the loop order of most pertur- 
bative calculations using lattice actions. In addition, we 
perform Monte Carlo simulations over a range of gauge 
couplings extending from the high-/? regime down to the 
lattice spacings used in typical simulations. We demon- 
strate that the two loop formula predicts the numerically 
self-consistent u^'^^) to within a few digits of the fourth 
decimal place and the additional tuning required is min- 
imal (especially when the action can be rescaled as in 
[4, 5]). For this reason we refer to both u^^"^^ and u^^*-^^ 
as u from hereonin. The small deviations at physical 
couplings are shown to be consistent with a third order 
correction to u and we infer the coefficient of this. In 
most cases no tuning at all is required if the third order 
is included. 

These calculations are carried out for two SU{?)) lat- 
tice gauge actions; the Wilson action and a first order 
Symanzik improved action. Isotropic and anisotropic 
lattices are studied and interpolations of the coefficients 
with the anisotropy are given. 

The structure of this paper is as follows. We first dis- 
cuss the perturbative calculation in Section II. In Sec- 
tion III we describe high-/? Monte Carlo simulations, and 
use the results to obtain higher order coefficients in the 
perturbative expansion, before concluding in Section IV. 
A brief description of twisted boundary conditions is 
given in Appendix A. 

The results presented here extend and, in some cases, 
correct the preliminary results presented in [6]. Exten- 
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FIG. 1: Feynman diagrams for the mean link. 
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FIG. 2: Feynman diagrams for the gluon self energy, E^„. 
Feynman rules are discussed in [5]. 

sion of this work to actions including fermions is being 
carried out and will be reported in a future publication. 



TABLE I: One loop coefficients for various anisotropies, x^ 
extrapolated to infinite lattice size, L. 



action 


anisotropy 






Wilson 


1.0 


-0.077467 


-0.077467 




1.25 


-0.087488 


-0.051452 




1.5 


-0.095096 


-0.036299 




2.0 


-0.105351 


-0.020351 




2.5 


-0.111543 


-0.012886 




3.0 


-0.115452 


-0.008851 




4.0 


-0.119871 


-0.004899 




6.0 


-0.123393 


-0.002142 




8.0 


-0.124713 


-0.001197 


Symanzik 


1.0 


-0.063049 


-0.064664 


improved 


1.25 


-0.071472 


-0.042877 




1.5 


-0.077940 


-0.030150 




2.0 


-0.086772 


-0.016934 




2.42 


-0.091467 


-0.011462 




2.5 


-0.092171 


-0.010720 




2.78 


-0.094281 


-0.008616 




3.0 


-0.095620 


-0.007365 




3.5 


-0.097922 


-0.005363 




4.0 


-0.099521 


-0.004077 




5.0 


-0.10152 


-0.002584 




6.0 


-0.10266 


-0.001784 




8.0 


-0.10384 


-0.000997 



II. PERTURBATION THEORY 



We write the Landau gauge mean links 

«^ = (trC/^) (1) 
(where N is the number of colours and tr = -^^ Tr) as 

n^ = l + a(V + a(V + C)(5')- (2) 

Using an anti-Hermitian gauge potential, = expg'A^, 
the mean Landau link in perturbation theory is 



(trC/^) = l + 5^/tr^ 



ir^\+0{g'). (3) 



We use twisted boundary conditions to regulate the gluon 
zero mode in a gauge invariant manner. A brief resume 
of relevant results and definitions is given in Appendix A 
and the reader is referred to [5, 7] for a fuller discussion. 



A. one loop 



The one loop contribution comes from 



1 1 



TW 



k,p 



(4) 



The one loop coefficient is thus the sum over the tree 
level propagator in Landau gauge (see Eqn. (18) of [5]), 



^ 2Vtw 



(5) 



shown in Fig. 1. We note here that the phase in the 
integrand is precisely that of the inverse propagator in 
the twisted formalism (see Eqn. (63) of [5]), so o^^^ is 
real. 



B. two loops 

The two loop contribution arises from two sources. The 
first is from the above integral, Eqn. (5), where the tree 
level propagator is replaced by the one loop self energy 
bubble shown in Fig. 2, giving a term 



a?) = 



2^1 



TW 



^^('='*=)G^,(fc)E,^(fc)G^^(fe) 



+ Eqn. (8). 



(6) 



As has the phase of the inverse propagator, this 
expression is also explicitly real. 
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TABLE II: Two loop coefficients for various anisotropies, Xi extrapolated to infinite lattice size, L. 



action 


anisotropy 










Wilson 


1.0 


-0.021093 (63) 


-0.021093 (63) 


0.002912 (63) 


0.002912 (63) 




1.25 


-0.023891 (74) 


-0.013222 (91) 


0.003573 (74) 


0.002930 (91) 




1.5 


-0.025721 (95) 


-0.008913 (63) 


0.004861 (95) 


0.002760 (63) 




2.0 


-0.027669 (91) 


-0.004730 (37) 


0.007771 (91) 


0.002116 (37) 




2.5 


-0.02827 (11) 


-0.002900 (23) 


0.01049 (11) 


0.001578 (23) 




3.0 


-0.02839 (12) 


-0.001916 (16) 


0.01262 (12) 


0.001228 (16) 




4.0 


-0.02811 (13) 


-0.001008 (14) 


0.01558 (13) 


0.000778 (14) 




6.0 


-0.02718 (27) 


-0.0004080 (42) 


0.01876 (27) 


0.0000377 (42) 


Symanzik 


1.0 


-0.01327 (23) 


-0.01435 (24) 


0.00273 (23) 


0.00206 (24) 


improved 


2.0 


-0.01831 (34) 


-0.003444 (63) 


0.00574 (34) 


0.001251 (63) 




3.0 


-0.01931 (50) 


-0.001456 (30) 


0.00882 (50) 


0.000711 (30) 




4.0 


-0.02011 (53) 


-0.0010482 (181) 


0.01001 (53) 


0.0001857 (181) 



The second term arises from the tree level portion of Ssi{Po,Xo,Us,Ut) 



the At term: 



1 1 



(tr^*) = T7T7^ Yl exp[i(fci + ... + A;4)-x] 

X Tr [r(fci)...r(fc4)] 4(fci)...4(fc4). (7) 
Using Eqn. (A9), this reduces to 
«<?> = Eqn. (6) 



TW 



k.p 



(8) 



Since the integrand is symmetric in k and p, only the real, 
cosine part of z^'^'P^ survives and the integral is real. 

As {k,p) depends only on the twist components of k, p, 
we note that if 



(9) 



is stored for each possible value of the associated twist 
vector (see Appendix for the definition of fe), the second 
expression in Eqn. (8) can be reduced to the product of 
two one-loop calculations. 



C. Numerical integration 

We consider the lattice Wilson action (W) and the 
Symanzik improved action (SI) defined in [8], both with 
tadpole improvement in the spatial (s) and temporal (t) 
directions. 



Sw{l3o,Xo,Us,Ut) 



Po Ps,s' 

Xo ^ , ui 



uiuf 



x,s>s' ^ * 



12 ul 

1 Rsa,t 



1 Rs's',. 

12 m6 



I 12 ; ' 



(10) 



where s, s' run over spatial links; Pg^s' and Ps^t are 1x1 
plaquettes; Rss,s' and Rss,t are 2x1 loops; and xo is the 
bare anisotropy. 

It is clearly inconvenient to work perturbatively with 
actions whose parameters are implicit functions of the 
expansion parameter g^. To circumvent this problem we 
define 

P^Po/{ulut), = glulut , X = XoUs/ut, (11) 
and find 

SwiPo, Xo, Us, Ut) = Sw{P, X, Us = Ut = I) , 
Ssi{Po,Xo,Us,Ut) = Ssi{/3,x,Us = Ut = I) + 
g^ASsi + 0{g^) , 
AS's/= rfsA5ct 

1 / [-Rss.s' + R. 



X 



+ xRss,t , (12) 



where dg is the one loop coefficient in the (self-consistent) 
perturbative expansion of the tadpole improvement fac- 
tor. For mean Landau link improvement ds = ai^'*. 

Perturbative expansions will initially be calculated as 
power series in g"^ and the expansion then re-expressed 
as a series in gQ. 

The Fcynman rules are obtained by perturbatively ex- 
panding the lattice actions using an automated computer 
code written in Python [6, 7]. The additional vertices 
associated with the ghost fields and the measure for Lan- 
dau gauge are given in [5, 7]. The Feynman diagrams for 
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(2) 

TABLE III: Contributions to the two loop coefficients, a I . The notation used in labeUing 


the contributions 


is explained in 


the text. 












action coefficient 


anisotropy 


VEGAS 


Figure of 8 


aV'ASct 


ASct 


(2) 

Wilson Os ' 


1.0 


-0.022969 (63) 


0.001875 








1.25 


-0.026283 (74) 


0.002392 








1.5 


-0.028547 (95) 


0.002826 








2.0 


-0.031138 (91) 


0.003468 








2.5 


-0.03216 (11) 


0.003888 








3.0 


-0.03255 (12) 


0.004166 








4.0 


-0.03260 (13) 


0.04488 








6.0 


-0.03194 (27) 


0.004759 








1.0 


-0.022969 (63) 


0.001875 








1.25 


-0.014049 (91) 


0.000827 








1.5 


-0.009323 (63) 


0.000410 








2.0 


-0.004859 (37) 


0.000130 








2.5 


-0.002952 (23) 


0.000052 








3.0 


-0.001941 (16) 


0.000025 








4.0 


-0.001015 (14) 


0.000008 








6.0 


-0.0004095 (42) 


0.000002 






(2) 

bymanzik Os 


1.0 


-0.01232 (23) 


0.0012422 


-0.0021920 


0.034766 


improved 


2.0 


-0.01693 (34) 


0.0023529 


-0.0037334 


0.043026 




3.0 


-0.01780 (50) 


0.0028573 


-0.0043696 


0.045697 




4.0 


-0.01855 (53) 


0.0030951 


-0.0046588 


0.046811 


(2) 


1.0 


-0.01468 (24) 


0.00130671 


-0.00097386 


0.015446 




2.0 


-0.003193 (63) 


0.00008961 


-0.00034017 


0.0039203 




3.0 


-0.00131 (30) 


0.00001695 


-0.00016325 


0.0017073 




4.0 


-0.0009590 (181) 


0.00000520 


-0.00009436 


0.0009481 



the one and two loop calculations are shown in Figs. 1 
and 2. 

The one loop integration is carried out by direct sum- 
mation of all twisted momentum modes for hypercubic 
lattices with 4^ < < 32. To speed up the approach 
to infinite volume, the momenta are "squashed" in the 
directions with periodic boundary conditions using the 
change of variables k ^ k' suggested by Liischer and 
Weisz [7] 

K^^i"" "a* sin(fc^) , (13) 

giving an integrand with much broader peaks which is 
easier to evaluate numerically. It is easy to see that a 
reasonable choice of parameter is ~ 1 — (xL^)~^ and 
significantly reduced the dependence on L. The calcula- 
tions were all possible on a workstation. 

All results were extrapolated to infinite volume using 
a fit function of the form cq -I- Ci/L^, which worked ex- 
tremely well for L > 16. As a further check, the coeffi- 
cients in the expansion of the spatial mean link should 
extrapolate to the same value in twisted and periodic di- 
rections. We found this to be the case, with the discrep- 
ancy much smaller than the number of significant figures 



quoted in the results in this paper. 

We calculate the two loop coefficient in three parts. 
The "figure of 8" diagram from Eqn. (8) is calculated by 
mode summation using the reduction described above. 

The remaining two loop calculations arising from 
Fig. 2{a — d) (and from the trivial Fig. 2e) were car- 
ried out using Vegas, a Monte Carlo estimation pro- 
gram [9, 10]. With a version of Vegas parallelised using 
the MPI protocol, these computations were carried out 
on between 64 and 256 processors of a Hitachi SR2201 
supercomputer. Each run took of the order of 24 hours. 
We label the results from this calculation "VEGAS". We 
used the 2 twist boundary conditions described in the 
Appendix. The lattice was of infinite extent in the un- 
twisted = 3, 4 directions, and of size L e {4, 8, 16, 32} 
in the twisted /i = 1,2 directions. Again, the Liischer- 
Weisz squashing was applied to the untwisted momentum 
components. 

Finally, for the Symanzik improved actions there is a 
counterterm insertion. Fig. 2/, arising from AS'ct, which 
requires a one loop calculation carried out by mode sum- 
mation. Multiplying this result by the appropriate dg 
gives the full contribution of the counterterm. 
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We show the one loop results in Table I and the two 
loop results for ds = a^s^^ in Table II. As it is possible that 
results for the mean link in Landau gauge from an action 
with a different tadpole improvement scheme might be 
desired (or an action with no improvement at all, when 
ds — 0), the contributions of the various parts of the two 
loop calculation are shown in Table III to allow the reader 
to construct the appropriate two loop contribution (the 
one loop numbers being unchanged) [4]. 



D. Resumming the series 

Perturbatively speaking, tadpole improvement resums 
the perturbative expansion in the hope of reducing trun- 
cation errors in the series. This amounts to using gQ as 
the expansion parameter, rather than g^. Writing 



ui = 



_jl + a^^g'+a^^g^+0{g') 
l + bWg2+b(^)gl+0{gl) 



(14) 



we find 



The coefficients b^^l for mean Landau link tadpole im- 
provement are given in Table II. 

Finally, we carry out interpolating fits in the 
anisotropy, x, to all perturbative coefficients using the 
function 



,«(3aW 



(15) 



ci 



C2 



Co H h ^ 

X X 



C3- 



logx 



(16) 



These fits work well. The fit parameters are given in 
Table III and are shown in Figs. 3 and 4. It should be 
noted when using these results that x is the anisotropy 
after the majority of tadpole factors have been scaled out 
of the action, as is shown in Eqn. (12). 



III. MONTE CARLO SIMULATION 

By comparing the truncated two loop perturbative se- 
ries obtained above with results from high-/? Monte Carlo 
simulations, higher order coefhcients in the perturbative 
expansion may be inferred [11]. 

A. Implementing tlie Langevin update 

Field configurations were generated using a 2nd or- 
der Runge-Kutta Langevin updating routine. The im- 
plementation of the Langevin evolution is such that any 
pure gauge action can be simulated by simply specifying 
the list of paths and the associated couplings. The group 
derivative for each loop is then computed automatically 
by an iterative procedure which moves around the loop 
link by link constructing the appropriate traceless anti- 
Hermitian contribution to the Langevin velocity. This is 
the most efficient implementation, minimising the num- 
ber of group multiplications needed and can be applied 
whenever the quantity to be differentiated is specified as 
a (closed) Wilson path. 

The method is as follows. We specify a given link by 
a start site, x, and a signed direction, ^: 



(17) 



0.02 



0.01 - 



-0.01 



-0.02 



-0.03 




0.02 



0.01 



-0.01 - 



-0.02 



-0.03 




J , L 

3 4 

X 



FIG. 3: Two loop coefficients to the Landau mean link for the 
Wilson action. 



FIG. 4: Two loop coefficients to the Landau mean link for the 
tadpole improved Symanzik action. 
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TABLE IV: Fits to perturbative expansion coefficients as functions of x, as defined in the text. 



action 



quantity 



const. 



1/x 



log X/X 



Wilson 



,(1) 



,(1) - 



-0.120648 




0.103142 
-0.0168915 



-0.0174148 


0.0296102 





0.0198738 
0.00384188 
-0.041532 
0.0180304 



-0.0599588 
-0.0609128 



-0.0235448 
-0.024975 
0.0148529 
-0.0152916 



-0.0614259 
0.00881346 



-0.0410654 
-0.00118467 
-0.0134185 
-0.00847315 



Symanzik 
improved 



,(1) 



-0.10043 




0.0915721 
-0.012287 



-0.0300447 



0.00122385 




-0.0215242 

0.00708635 
-0.0728384 
0.0129186 



-0.0541941 
-0.0526966 



0.038299 
-0.0214381 
0.0743446 
-0.0108602 



-0.0536973 
0.00656986 



0.0372853 
-0.00419861 

0.0644862 
-0.00676255 



The for ^ > are vectors of the unit cell and = 
—e-n for negative /z. 

The Wilson path, W{x,U), of length M is a unitary 
matrix specified by a starting point x (where it is ini- 
tially stored) and an ordered list of M links in directions 
(/Ui, ...,/iM). Explicitly, 

Wix,U)= II U^^ix + h-i), (18) 
with /m = (1, 2, . . . , M) and 

r 0, 2 = 0, 

Ee..> 0<^<M-l, (19) 

The contribution to the pure glue action from closed 
loops is then 

S{W) =gwY.Re Tt{W{x, U)) , (20) 

X 

where gw is the associated coupling constant and 
W{x,U) contributes to the Langevin velocity v{xj,U) 
at all sites xj = x + Ij for j = ... M — 1 

v{xj, U) = gwA[Wj{xj,U)] , (21) 

where the operator A projects onto the SU{N) algebra 

of anti-hermitian, traceless matrices. If is the j-th 
left cyclic permutation of Im, 

Wj{xj,U)= H U^,{x + li-i). (22) 

Clearly, Wj{xj, U) can be obtained from Wj-i{xj-i,U) 
stored at site Xj^i by only two matrix multiplies and 
then the result is stored at site Xj. In principle, the 



whole lattice can be treated in parallel. The outcome 
is that the number of matrix multiplications required to 
compute V is 3(M— 1) rather than the M(M— 1) for naive 
construction. This improvement is particularly relevant 
in application to improved QCD actions with dynamical 
fermions. 

The 2nd order Runge-Kutta algorithm used is a two 
step update: 

J/«(r,a;) = R^^\T,dT,x)U,{T,x) , 
U^{T + dT,x) = R'-^\T,dT,x)Uf,{T,x) , 
R^^\T,dT,x) = ex.p{u^^\T,dT,x)) , 
R^'^\T,dT,x) = exp(u^^^(r, rfr, x)) , 
u'^^\T,dT,x) = dT/2v{x,U{T)) + 

Vdrrj'^^^TjX) , 
u'^^\T,dT,x) = dTC+v{x,U^^\T)) + 

y^d^ ir]^'\T,x)+r]^^\T,x)) . 

(23) 

Here, C± = 1 ± with Cat being the adjoint 

Casimir for SU{N) {C3 = 3). The r]^P^T,x) for p = 1,2 
are N x N anti-hermitian Gaussian random matrices: 

V^\t,x)=tj^J'\t,x)T,, (24) 

[with Ta being the (anti-hermitian) generators of 

SU{N)]. They satisfy 

( ) (r, X) V^f\T', X') ) = V Srr' ^xx' ^aa' ■ (25) 

B. Twisted BC and tunnelling 

All 2-, 3- and 4-twist boundary conditions were im- 
plemented for simulation. The manner of the implemen- 
tation is as described in [7]. An explicit representation 
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of the Q matrices is not needed for the update for the 
pure gauge (quenched) update since the field simulated 
is U which is related to the perturbative field U by the 
redefinition 



u^{x) = u^{x)n^{x) 

j 1 Xi^T^Li^ 



(26) 



The action in terms of the {[/} is identical to the un- 
twisted action except that a loop whose projection onto 
the (p, i/)-plane (p < v) encircles the point {L^ + 
1/2, Lu + 1/2) has an additional factor of (z^^i/)^"^ where 
c is the integer winding number [i.e. the sense of the 
circulation: c > (c < 0) for anti-clockwise (clockwise) 
circulation of |c| turns]. This weights the corresponding 
contribution to the action by the factor cos {2-Kcnfj,v/N) 
compared with the untwisted case. An explicit repre- 
sentation of ri would be needed for a simulation with 
dynamical fermions (and, as we shall see, for gauge fix- 
ing). 

For twisted boundary conditions the number of differ- 
ent Z{N) phases is N"^ [12] and these are distinguished 
by order parameters which arc the set of independent 
non-zero Polyakov lines winding around the lattice. For 
example, one such Polyakov line in the 4-twist case is 
-Pi32(a;) where 

i,; 

Pi{x,U) = Ui{x + TliBi) , 

ni=0 

Pijk{x,U) = Ti[Pi{x,U)Pj{x,U)Pk{x,U)] . (27) 

Here i,j,k are signed directions as above. Tunnelling 

can occur between different Zj^ vacua, corresponding to 
the miiltiplication of all links on a space-like slice by an 
element z from the centre group of SU{N). 

Twisted boundary conditions create a barrier between 
the different Z{N) phases and the zero mode, responsi- 
ble for non-perturbative effects [13], is absent. This is a 
prerequisite for comparison with perturbation theory and 
the measurement of coefficients of terms in the series of 
higher order than those calculated. The study in refer- 
ence [11] shows that with 3-twist boundary conditions 
tunnelling is absent for /3 > 9, even on small lattices of 
4*. The probability for tunnelling then reduces as the 
lattice volume increases. We use a 8'^ x 16 lattice and 
with 4-twist boundary conditions we expect tunnelling 
to be rare or absent for lower values of (3. We can obtain 
an idea of whether it is occurring by plotting the tempo- 
ral history of a non-zero Polyakov line such as P132. The 
results are shown in Fig. 5 for f3 = 9.0 and f3 = 6.0. The 
density of points in both cases is centred on a point on the 
real axis and no points are found in the centre-equivalent 
positions on the lines in directions (— 1/2, ±-\/3/2). We 
conclude that tunnelling events are absent and that we 
can be confident that the influence of non-perturbative ef- 
fects on our measurements has been minimised. For this 



reason, and on account of the convergence of the gauge 
fixing discussed below, we restricted our simulations to 
the 4-twist boundary conditions only. 



C. Fixing to Landau gauge 

To fix each configuration to Landau gauge we maximise 
the corresponding gauge function with respect to gauge 
transformations. In the continuum this function is 

J^{{A}) = J d^x tvA^{x)A^{x) , (28) 
and the lattice analogue that we use is 



, (29) 



where the c„ arc coefficients chosen so that terms in 
the perturbative expansion of J^l between and some 
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FIG. 5: The scatter plot for the expectation for Polyakov 
line Pi32 for /3 = 10 (top plot) and /3 = 6 (bottom plot). 
Although the mean value is becoming small as /3 decreases, 
both distributions are centred on a point on the real axis and 
show no sign of tunnelling to centre-equivalent distributions. 
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higher power are absent; and /i,2,3 = 1/X) fi = X- 
practice, we use 



Cl = 1, C2 



1 



1 

63 



, C4 



1 

896 ' 



(30) 



with c„>4 = 0. This ensures that the definition of Landau 
gauge agrees with that used in the perturbation theory up 
to and including 0{g^). This allows us to check the two 
loop perturbation theory for Ug, Ut and to deduce higher 
loop contributions using the simulation. Our improved 
gauge-fixing function is not the same as those described 
elsewhere [14]. It docs, however, correspond to the defi- 
nition of Landau gauge used in the perturbation theory, 
namely 



A(,-)^^(a;^ + e^/2) = 0, 



(31) 



where A^(x + e.n/2) is the pcrturbative gauge field. 

A Fourier accelerated algorithm is vital when fixing 
to Landau gauge, which is otherwise prohibitively time 
consuming. To carry out the Fourier transform on fields 
with twisted boundary condition we define r(Ti, a;) by 



r(n, x) = r(n) exp(27rm^a;^/A'') 



(32) 



where the constant matrices r(n) are defined in the Ap- 
pendix. For a twisted field v{x) in the algebra of SU{N) 
we define 



j{n,x) = Tr [r^{n,x)v{x)] 



(33) 



Now v{n,x) is periodic on the lattice and we can de- 
fine the periodic Fourier transform of v{x), which can be 
implemented by the FFT algorithm, by 

v{n, k) = v{n, x) exp(— 27rifc^a;^/7V) , (34) 

X 

where fc^ is given in terms of and in Eqn. (A5) 
and hence the argument n on the LHS is redundant and 
retained only as a reminder: n can be deduced from k 
but it is convenient to explicitly display the dependence 
of V on n. The inverse transform is then 

v{n,x) = — ^ v{n, k) exp(27rifc^a;^/iV) , 



TW 



v{x) = r(n, x)v{n, x) , 



(35) 



(see [5] for further details). 

To maximise J^l{{U}) with respect to gauge transfor- 
mations we define the gauge velocity 



i0{x) = Ta 



5ea{x) 



J-4{C/('^^})|.=o , 



(36) 



and use steepest ascent. The gauge transformation g{x) 
is given by g{x) = expeo(a;)T(j and satisfies the twisted 
boundary conditions 



g{x + L^e^) = nij,g{x)n'l, 



(37) 



and then 



Ujf\x) = gix)U^{x)g\x + e^) 



(38) 



We, however, simulate with the periodic field {U} which 
transforms under gauge transformations as 



m<'\x) = g{x)%{x)g'<{x + e^) 



with the rule 



g{x) = g{x), g(x + Lf,ef,) = g{x) 



(39) 



(40) 



Note that, whilst the {U} transform effectively with a 
periodic gauge transformation, the infinitesimal gauge 
transformation still has the twisted spectrum derived 
from Eqn. (37). The outcome is that to maximise 
J^l{{U}) we apply successive infinitesimal gauge trans- 
formations g{x) = exp{eu){x)) to the configuration {U} 
where 

u;{x) =(3Y,nCr,f^A[{U^{x)r - {Ul{x - e^))"] , 

n^fx 

(41) 

with Ufi{x) = Ufj_{x)Cljj^(x) and e is a small step size. The 

simulation uses the fields {U} and so an explicit repre- 
sentation of the ilf_i is needed (such as in the Appendix). 

The Fourier acceleration of gauge fixing has been dis- 
cussed in [15] and the Fourier accelerated gauge velocity 
is given by 



uiaccin, k) = K{k) u){n, k) 



(42) 



where the tilde notation is defined in Eqns. (33-35). The 
acceleration kernel, K{k), is given by 



K{k) 



fc2 



(43) 



where m is a mass parameter which may be set to zero 
for twisted momenta as there is no zero mode. Note 
that since k^ depends on n for given periodic component 
fe, then K{k) distinguishes between modes labelled by 
different n. This is significant only for small k where, 
nevertheless, the acceleration has the largest effect and 
careful control of the IR behaviour of the algorithm is 
needed. 

The gauge transformation g{x) = exp(e[iacc(a^)) is 
constructed for small step size e. Because the expo- 
nentiation and Fourier acceleration are expensive it is 
best to successively apply this gauge transformation to 
the configuration {tl} until a maximum of J-" on the 
given gauge orbit is achieved and then cDacc is recom- 
puted and the procedure repeated. To speed up the pro- 
cess a sequence of gaiige transformations is constructed, 
g„i = .9m- 1 1 .91 = ffj m = 1, . . . , M. The maximum on 
the given orbit is determined using gM and then succes- 
sively refined using the gm in descending order in m. This 
reduces the number of matrix multiplies required needed 
whilst allowing e to be chosen to be very small. We chose 
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M = 5. The algorithm is applied until the condition on 
the Fourier-accelerated velocity 




is satisfied. We found that S ^ 10~^ to be more than 
sufficient to fix the to Landau gauge very accurately. 

We observe that (a) without Fourier acceleration con- 
vergence is prohibitively slow, and (b) that convergence 
was improved by at least a factor of 10, measured in com- 
puter time, when 4-twist boundary conditions were used 
compared with the 2- or 3-twist boundary conditions. 
The reason for the latter result is unclear. It is possi- 
ble that tunnelling between different toron vacua is the 
most strongly suppressed in the 4-twist case, but the im- 
provement was still much in evidence at large /?, /3 ~ 30, 
and from the mean Polyakov line in the untwisted di- 
rection(s) there was no noticeable signal that tunnelling 
was occurring in the 2- and 3-twist cases. The number 
of near zero modes is fewer in the 4~twist case but there 
is no convincing argument that the effect should be so 
dramatic. 



D. Gribov copies 

Fixing {[/} uniquely to Landau gauge corresponds to 
finding the global maximum of !Fl{{U}). Other (local) 
maxima give rise to lattice Gribov copies of this gauge. 
The global maxima lie in the Fundamental Modular Re- 
gion (FMR) of the Landau gauge cross-section and since 
the FMR contains the identity we expect that our per- 
turbative calculations relate to fields gauge-fixed in this 
way. The mean link values are closely related to J^l at 
the global maximum and to use Gribov copies to evalu- 
ate them will clearly give, on average, a lower value for 
the tadpole coefficients than predicted by the perturba- 
tion theory. We have not investigated the role of Gribov 
copies in the simulation and so, in principle, the mea- 
sured tadpole coefficients will have a negative systematic 
error. However, we believe that the effect is negligible for 
large enough /? and when fully twisted boundary condi- 
tions are employed. We defer further discussion of this 
point until the next section. 

E. Simulations 

The simulations were carried out for a range of (3 values 
from 5.2 < /3 < 30 for x = 2 which easily encompasses 
the physical region. Because it is more physically rele- 
vant we concentrated on the improved Symanzik action. 
The simulations were run on 64 processors on the Hitachi 
SR2201 computer at the Tokyo Computer Centre on lat- 
tice size 8^, and on 48 processors on the SunFire F15K 
of the Cambridge-Cranfield High Performance Comput- 
ing Facility on an 8^ x 16 lattice. As mentioned above. 



we used the 4-twist boundary conditions in these calcu- 
lations. 

Using the actions specified on the RHS of Eqn. (12) 
the self-consistent tuning of {us. ui) is fast. In fact, it is 
immediate for the Wilson action since there is no explicit 
reference to {us,ut) once the rescaling in Eqn. (11) has 
been done. For the Symanzik action there is a residual 
dependence through the counter-term ASsi on Us but 
starting at 0{g^). In this case, in the simulation a value 
for Us is chosen and the mean spatial link is measured and 
the next value of Ug inferred which is used as the starting 
value for a new simulation. This process needs to be 
iterated only two or three times for accurate convergence 
to the self-consistent values of {us,ut) to be determined. 
For more on the efficient self-consistent timing of tadpole 
improvement factors in Monte Carlo simulations, see [4] 
and the linear map techniques in [16]. It is, of course, an 
aim of this study to remove the need for such tuning. 

Precision comparison of perturbative and simulation 
mean links requires consideration of finite size effects 
(FSEs), and we would like to compare the two for a given 
volume. From looking at the variation with L of the per- 
turbative results, the major finite size effect was found 
to be from the one loop, O^g^) contribution and the sim- 
ulation agreed accurately with the L = 8 prediction. Al- 
ready for L > 8 the higher-order coefficients show very 
little sensitivity to L and so for comparing with simula- 
tion we set them at their large-L values. We then com- 
pare the results of the simulation with the calculated two 
loop expansion. In Fig. 6 we plot 

u_=«(^^)-(l + 6[^g + 6p)5o^) , (45) 

the measured mean link with the two loop perturbative 
prediction subtracted, against for x = 2. 

The object is to show that the simulation is very ac- 
curate and agrees very well with the perturbative pre- 
diction for small enough g^ and to also infer the size of 
the deviation at larger gg due to three loop and higher 
order terms, i.e. at and higher. A polynomial fit will 
give a good estimate for higher coefficients in the self- 
consistent perturbation series. Care must be taken when 
fitting polynomials to data in a given range of the expan- 
sion parameter. If the true coefficients are such that in 
this finite window there is an approximate cancellation 
between some combination of terms, then all such terms 
must be included in the fit function. If not, as might 
be the case when trying polynomials of increasing order, 
marked instability in goodness of fit and even low order 
coefficients will be seen. Inclusion of all orders can be 
approached using 'Bayesian' techniques [17]. These did 
not, however, appear necessary in describing our data, 
and the results quoted come from straight polynomial 
fits at the given order. 

Despite comparing simulations with the L = 8 pertur- 
bation theory, there can, however, still be slight FSEs as 
tunnelling has forced us to use different boundary condi- 
tions in the simulation and perturbation theory (which 
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matters at finite L). To allow for this small discrepancy, 
we fit polynomials of the form 

u-=5b['^9^ + bf'^gl, (46) 

where Sbl ' is the finite size correction at two loops. The 
quality of the fits is good down even to gauge couplings 
/? 5 — 6, with no sign of terms of C(5o)- The x^/d.o.f 
are both 0.6 and the fit coefficients are 

Sb'^^^ = 0.00032(4) , 6(3) = -0.00097(3) , 

56f ) = 0.00015(1) , 6f ) = -0.00029(1) . (47) 

The finite size correction corresponds to a 2% eflfect in 
the two loop calculation. This is certainly very reason- 
able. The measurement of the three loop coefficients is 
accurate to 3% and we should also expect a finite size 
error of a similar order. 



IV. DISCUSSION AND CONCLUSIONS 

Tadpole improvement via Landau mean link factors is 
becoming increasingly important in lattice field theory, 
but there are computational costs in tuning the mean- 
link factors to their self-consistent values. In some cases, 
such as the small pure gauge theory lattices used in this 
work, the costs are bearable. In other cases they become 
significant. For example, highly anisotropic lattices, such 
as would be used for finite temperature studies, would 
need to be tuned at zero temperature for a physical spa- 
tial volume. The very large time extent then required 
makes numerical tuning a major study in itself. 

In this paper we calculate the mean link in Landau 
gauge perturbatively to two loop order. By comparing 

0.002 I ^ ; ^ ; ^ ; ^ , 

C 
O 



CO 




So 



FIG. 6: Simulation of tlic Landau mean link tadpole parame- 
ters for the Symanzik-improved action with x = 2 versus self- 
consistent coupling with second order perturbation theory 
subtracted. The solid lines show fits with three loop, O{go) 
parametrisation. For fit results see text. 



this to simulations, we have shown that this result pre- 
dicts the mean link in Landau gauge to a high accuracy. 
Use of the formulae obtained remove the need for the 
expensive tuning. 

We have used the automatic generation of the pcrtur- 
bative vertices to calculate the self-consistent tadpole- 
improvement parameters Us, ut for the anisotropic Wil- 
son and Symanzik-improved gluon action to two loops as 
a function of the lattice coupling and the anisotropy 
coefficient x- The tadpole-improved actions, Eqn. (10) 
are parametrised by go and xo = XUs/uf The one 
and two loop coefficients are given in Tables I and II, 
respectively, for different %• The calculation was done 
in Landau gauge with twisted boundary conditions in 
the /i = 1,2 directions so that there was no zero mode. 
Whilst the one loop and figure-of-eight integrals were 
done with mode summation on an lattice, the ma- 
jor two loop numerical Feynman integrations were done 
using parallel Vegas [9, 10] on an x oo^ lattice. For 
L > IQ (with appropriate "squashing" of the momentum 
variables [7]) there was no observable i-dependence. The 
coefficients of and 5* for the direct perturbation series 
and of 90 and for the self-consistent perturbation series 
are given in Tables I-III for different anisotropics, x- De- 
tails of the component calculations are also given in these 
tables to help facilitate the future calculation of different 
but related quantities. Interpolating fits were obtained 
to these coefficients as a function of x (Table IV). 

The simulation was done using a second order 
Langevin algorithm for the field update and the Lan- 
dau gauge is fixed for each configuration using a Fourier 
accelerated steepest ascent method to maximise the (im- 
proved) gauge function J-{U). 

We found that fixing to Landau gauge was very slow 
unless Fourier acceleration was used, and that by using 
twisted boundary conditions in all four directions (4- 
twist) the convergence of the algorithm was improved 
by at least a factor of ten compared with applying the 
twist in just two or three directions (2- or 3-twist). The 
reason for this is not clear but may have a connection 
with the suppression of tunnelling between different toron 
vacua, although no clear signal supporting this surmise 
was foimd. Practically, this effect was welcome in the 
isotropic case, and essential for x = 2 where the conver- 
gence was otherwise prohibitively slow. 

We did not investigate the role of lattice Gribov copies 
in fixing to Landau gauge by maximizing !Fl{{U}) in 
Eqn. (28). There has been a long discussion in the lit- 
erature concerning the role of Gribov copies [18-23]. In 
principle, the effect of selecting a local rather than the 
global maximum of will be to reduce the measure val- 
ues of the tadpole improvement coefficients. We believe 
that the effect in this calculation is negligible for a num- 
ber of reasons. We do not expect a problem with Gribov 
copies at large-/? but some effect may occur as the physi- 
cal region is approached [24] . However, the use of twisted 
boundary conditions is likely to suppress the discrepancy 
due to the Gribov problem. This is because there is no 
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zero-mode with such boundary conditions and by using a 
twist in aU four directions the number of low-lying modes 
is reduced. In [20] the effect of the zero-mode was shown 
to be very important in the U{1) case and by performing 
the extra gauge transformation to remove its effect, the 
ZML gauge, the distribution of mean link values mea- 
sured was greatly reduced. Another effect is due to dou- 
ble Dirac sheets (DDS) which have zero action but are 
strongly correlated in the U{1) theory with a reduction 
in the mean link value and which are mooted to also be 
relevant in SU (N) theories. We see no evidence for DDS 
contributions in that we see no negative spikes of the na- 
ture clearly observed in the U{1) case [20]. However, we 
have no quantitative analysis to support this remark. 

If there are Gribov effects in our simulation then these 
will be included in the spread of results as a systematic 
error and are therefore to a great extent included in our 
quoted error which we perforce analysed by purely sta- 
tistical methods [25]. 

The Landau gauge function is improved; it is cho- 
sen so that it corresponds to the lattice Landau gauge 
used in the perturbative calculation up to term of or- 
der 90°. Thus, perturbation theory and the simulation 
should agree up to and including four loops. 

The simulations were done on an x T lattice with 
L = 8, T = 16 for the Symanzik-improved action with 
X = 2 whch corresponds to an equal-side hypercubical 
lattice when measured in physical units. The couplings 
used were in the range 5.2 < (3 < 30. The simulation 
boundary conditions were twisted in all four directions, 
and in order to compare with calculation an adjustment 
for finite L,T effects is needed. Finite size effects were 
only appreciable in the one loop coefficients. Accordingly, 
we used one loop coefficients calculated on a finite lattice, 
and two loop data for infinite volume. The higher or- 
der perturbative coefficients that are extracted may thus 
change slightly for larger L. 

We fitted the data (less the two loop perturbative re- 
sult) with a polynomial containing terms in and g^. 
This described the data extremely well over the full range 
of gauge couplings. The first term in the fit function al- 
lowed for a small finite size correction to the two loop 
perturbative result. This finite size discrepancy corre- 



sponds to a 2% correction to the two loop calculation. 
This is certainly very reasonable. The measurement of 
the three loop coefficients is accurate to 3% and we should 
also expect a finite size error of a similar order. The main 
conclusion is that even in the physical region /3 ~ 5 — 6 
perturbation theory works very well indeed for the Lan- 
dau mean link and, moreover, there is no observable de- 
viation from a three loop, C(,9o) perturbative approxi- 
mation. This is very encouraging for the accurate design 
of QCD actions and the corresponding perturbative anal- 
yses based upon them. 

In principle, a full simulation analysis would give three 
loop numbers for a range of x values, enabling a fit of 
the same form as is given in Eqn. (16). However, these 
coefficients are seen from the figine to be small even in 
the physical region, and so we conclude that Ug and Ut 
are approximated sufficiently accurately by the two loop 
expressions for all practical purposes. 

Finally, we remark that it could be possible to calculate 
Us and Ut by simulation using stochastic perturbation 

theory. However, although it has proved very successful 
for the expansion of the mean plaquette, the gauge fix- 
ing required here makes this approach extremely slow. 
Nonetheless, in general, the double expansion of the ac- 
tion and then the gauge potential itself as power series in 
the coupling required for the stochastic evolution equa- 
tions is a task to which the Python code can be readily 
adapted. 
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APPENDIX A: TWISTED BOUNDARY 
CONDITIONS 

We give here a brief description of the boundary con- 
ditions. Further details can be found in [5, 7]. 

For an orthogonal twist the twisted boundary condi- 
tion for gauge fields (and potentials) is 



(Al) 



where the twist matrices fl^, are constant SU {N) matri- 
ces which satisfy 



(A2) 



and Zf^i, = exp{2'Kinij,v/N) is an element of the centre 
of SU {N) . The particular boimdary conditions imposed 
are uniquely specified by the antisymmetric integer ten- 
sor n^^, n^u e -E]v [12]. The twisted boundary con- 
ditions can be chosen to apply in two directions only, 
taken here to be the (1,2) directions, but can also be 
applied to the 3 and 4 directions if required. Our imple- 
mentation of these boundary conditions correspond to 
orthogonal twists since k = e^^apn^iyUap = mod N, 
and consequently only configurations with integral topo- 
logical charge, Q, will occur [12]; the perturbation theory 
is then correctly associated with the Q = sector. In this 
case fis, ^4 can be expressed in terms of fii, 0,2 once the 
values of n^^ are given. We use 



(A3) 
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The different choices of boundary condition are then 
given by assignments to fl^ and 0.4 as follows: 



0.4 (r^l2,r^l3,nl4,n23,n24,T^34) 



1 1 

o\ol 1 
o\ol OlOl 



(1,1,0,0,0,0) 
(1,1,0,1,0,0) 

(1,1,1,1,1,1) 



(A4) 



We refer to these choices as 2-,3-,4-twist boundary con- 
ditions, respectively. 

If the lattice has extent in the /U-direction, the mo- 
mentum spectrum is 



27rfc^ 27rn^ 
0<n^<N , < fc^ < . 



(A5) 



For 2-twist boundary conditions, 713 = 714 = 0. The 
3-twist boundary conditions are given by — ~{ni + 
TI2), 714 = and the 4 -twist ones by ns — —{ni + ^2), 
724 = rii — 712- In all cases we exclude the zero mode 
Til = 7l2 = 0. 

Negative momentum in these directions is — fc^ = 
adding appropriate multiples of L/j^ and 
respectively to remain in the ranges defined above. 

Twisted boundary conditions thus replace the N"^ — 1 
colour components of an SU (N) field by a similar number 
of extra momentum components, interstitial to the usual 
reciprocal lattice. The Fourier expansion of a gauge po- 
tential is now 

M^) = E e"--r{n)A^{k) , (A6) 

where Vtw = ^ is the twisted volume of the lat- 

tice. The sum over momentum indicates a sum over all 
and the components of the twist vector n = (tii, 712)- 
The — 1 SU{N) twist matrices are given in terms of 
Oi^2 by 

r{n) = ^("i+"2)("i-Fn2-l)/2Q-n2Qni ^ ^^j^ 

where z = ex]){2i'K/N) is an element of the centre of 
SU{N). 

We need only the trace algebra associated with these. 
Defining the symmetric and antisymmetric products of 
twist vectors 

(n,m) = 7li772i -I- 7l27n2 -I- (7I1 -I- 772i)(7l2 -I- 7712) , 
(n, m) = 7117722 — 7127711 , (A8) 

we obtain 

T{n) = 1 n^O mod N , 

Tr [r(n)] n^O mod A'' , 

r(n)t = z-5("'")r(-n) , 

r(n')r(n) = z5«"''">-("''"))r(n'-hn) , 

Tr [r(n')r(n)] = Ar05(-'-),5„,_„, , (A9) 
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where (n mod N) is understood to apply to each com- 
ponent, ni.2, as is the (5-fimction. The argument of F is 
evaluated mod N. 

Note that for clarity in the main text we some- 
times replace the twist vectors in expressions such as 
Eqns. (A8,A9) with their associated momenta. In each 



case, the argument is understood to be the twist vector 

alone. 

Note also that in Eqn. (A6) the location of the gauge 
field, X, is taken to be midway along the associated link. 
The components, often then expressed as integer 

multiples of units of half a lattice spacing. 



